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Abstract 


The implementation of the k-kL turbulence model using multiple computational fluid dy- 
namics (CFD) codes is reported herein. The k-kL model is a two-equation turbulence model 
based on Abdol-Hamid’s closure and Menter’s modification to Rotta’s two-equation model. 
Rotta shows that a reliable transport equation can be formed from the turbulent length scale 
L, and the turbulent kinetic energy k. Rotta’s equation is well suited for term-by-term mod- 
eling and displays useful features compared to other two-equation models. An important 
difference is that this formulation leads to the inclusion of higher-order velocity derivatives 
in the source terms of the scale equations. This can enhance the ability of the Reynolds- 
averaged Navier-Stokes (RANS) solvers to simulate unsteady flows. The present report 
documents the formulation of the model as implemented in the CFD codes Fun3D and 
CFL3D. Methodology, verification and validation examples are shown. Attached and sepa- 
rated flow cases are documented and compared with experimental data. The results show 
generally very good comparisons with canonical and experimental data, as well as matching 
results code-to-code. The results from this formulation are similar or better than results 
using the SST turbulence model. 


Nomenclature 

Roman letters 

Cd drag coefficient 

% Heaviside function 

L turbulent length scale 

L v k von Karman length scale 

M Mach number 

M t turbulent Mach number, y^k/a 2 

N number of nodes in grid 

Pk production of turbulent kinetic energy 

Re Reynolds number 

Re L Reynolds number based on length L 

Reg Reynolds number based on momentum thickness 

S, S ij symmetric velocity gradient tensor 

U,itj Cartesian velocity vector, ( u,v,w) T 

a local speed of sound 

b half width of stream, where (u — u e ) is half of (u m — u e ) 

c f local skin friction coefficient 

d distance normal to surface 

f c auxiliary function in compressibility model 

f$ auxiliary function in kL transport equation 

h grid spacing measure 

k turbulent kinetic energy 

p pressure 

r radius 

T temperature 

t time 
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u e velocity at edge of outer planar shear stream 
u m peak velocity of planar shear stream 
u + velocity in wall units 

Xi Cartesian coordinates, (x, y, z) 

y + distance normal from surface in wall units 

Subscripts 

acoustic based on ambient conditions 
oo free-stream condition 

jet relating to jet condition 

max maximum 

min minimum 

splitter splitter plate in planar shear case 
t total condition 

Conventions 

2D two dimensional 

ARN acoustic research nozzle 

CFD computational Fluid Dynamics 

C-D convergent-divergent 

EASM explicit algebraic Reynolds stress 

LES large eddy simulation 

NPR nozzle pressure ratio, Pt.jet/Poo 

NTR nozzle temperature ratio, T t j et /T 0 0 

RANS Reynolds averaged Navier-Stokes 

SAS scale-adaptive simulation 

SST shear stress transport 

URANS unsteady Reynolds averaged Navier-Stokes 
Symbols 

Sij Kronecker delta 

e scalar dissipation 

9 momentum thickness 

k von Karman constant 

fi bulk viscosity 

Pt turbulent eddy viscosity 

ui specific dissipation rate 

$ kL 

p density 

Tij Reynolds stress tensor 

Superscripts 

' first derivative 

" second derivative 
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1 Introduction 


W HILE two-equation models have been used routinely for the last 50 years, they are 
based on the Reynolds stress equations and, typically, a modeled length scale. The 
mechanism of the second equation for determining turbulent length scale is not fully un- 
derstood and a number of formulations use a special boundary condition for simulating 
its wall boundary condition. Even the more complex model closures like Reynolds stress 
models (RSM), or explicit algebraic Reynolds stress models (EASM) still use a length scale 
equation based on an underlying two-equation model. Almost all two-equation models use 
the turbulent kinetic energy, k, and its transport equation as one of the primary variables. 

Historically, the modeling of the second equation using dimensional arguments has been 
purely heuristic [1]. Many of the models for the production are only using strain-rate or 
vorticity derived from the mean flow terms, resulting in only one scale from the equilibrium 
of source terms for both equations. The second equation is considered, in most cases, 
the weakest link in turbulence models, including much more complex approaches such as 
differential Reynolds stress and hybrid RANS/LES formulations. It is difficult to justify 
using any of the complex turbulence models without fixing or using an alternate form for the 
second transport equation. One of the few exceptions is the modeling concept proposed by 
Rotta [2], which can be formed as an exact transport equation for the turbulent length scale, 
L. Rotta’s approach is well suited for term-by-term modeling and displays very favorable 
characteristics, as compared to other approaches. A key difference is the inclusion of higher- 
order velocity derivatives in the source terms of the scale equation. This potentially allows 
for resolution of the turbulent spectrum in unsteady flows. 

Menter et al. [3-5] presented a complete form of the k-\/kL two-equation turbulence 
model based on the Rotta [2] approach. In Menter, it was proposed to replace the problem- 
atic third derivative of the velocity, that occurred in Rotta’s original model, with second 
derivatives of the velocity. Menter utilized this two-equation turbulence model to formulate 
the Scale-Adaptive Simulation (SAS) term that can be added to other two-equation mod- 
els, such as Menter’s SST [6]. The SAS concept is based on the introduction of the von 
Karrnan length scale into the turbulence scale equation. The information provided by the 
von Karrnan length scale allows SAS models to dynamically adjust to resolved structures 
in unsteady RANS (URANS) simulations. This can create LES-like behavior in unsteady 
regions of flowfields. At the same time, the model provides standard RANS capabilities in 
stable flow regions. 

Abdol-Hamid [7] documented an initial form of the k-kL two-equation turbulence model. 
The reference showed the process to calibrate the constants within the range suggested by 
Rotta [2] and satisfying the near-wall logarithmic requirements. This model does not use a 
blending function to merge two turbulence scale equations as is done with the SST model, 
or have the near-wall damping functions typical of k-e turbulence models. It naturally 
contains the SAS characteristics through the von Karrnan length scale. The basic model 
was implemented in the CFD code PAB3D [8]. The formulations, usage methodology, and 
validation examples were compiled and presented to demonstrate the capabilities of the 
model. The model provides proper RANS performance in stable flow regions and allows the 
break-up of large turbulent structures for unstable flow regions, for example a cylinder in 
cross-flow or flow in a cavity. 

In the present report, we document the k-kL model and the process to implement it 
in CFD solvers. The model (with a few very minor differences from the model originally 
implemented in PAB3D) is in both the unstructured code Fun 3D and the multiblock struc- 
tured code CFL3D. This new model version is named k-kL-MEAH2015. This report shows 
the comparisons between the results generated from both codes to validate and verify the 
implementations. Results for all near-wall flows have been calculated using the grids that 
resolve the viscous sublayer with y + < 1. Simulations have been carried out using 5 grid 
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levels for the verification cases, avoiding grid refinement uncertainties. The validation re- 
sults are compared with available experimental and/or theoretical data, depending upon 
the case. Most of the cases are taken from the turbulence modeling resource website [9,10]. 

Subsonic and supersonic jet flows are quite difficult to predict with most R.ANS turbu- 
lence models. For subsonic jets, most turbulence models incorrectly predict the mixing rate 
so that the jet core length differs significantly from what is physically observed. The k-kL 
model also predicts a core length that is too short. A proposed jet correction, including 
compressibility effects, is described and is designated k-kL-MEAH2015+J. 


2 Turbulence Model Description 

The baseline k-kL turbulence is described in 2.1. Models to correct for free shear flows and 
compressibility effects are described in section 2.2. 


2.1 Baseline k-kL Model 


The k-kL-MEAH2015 two-equation turbulence model, Eqs. 1 through 12 , is based on the 
approach of Menter [3-5] to develop a k-\/kL model. A complete list of coefficients used 
by the present model is defined. The model is based on Rotta’s k-kL ($ = kL) with the 
modifications proposed by Menter where the third derivative of velocity was replaced with 
the second derivative of velocity. The closure constants were derived and documented by 
Abdol-Hamid [7]. The current version, k-kL-MEAH2015, is slightly different from the earlier 
version, k-kL-MEAH2013, that was reported by Abdol-Hamid [7]. 

d -> k 5 / 2 k 

—pk + V • (pUk) = P k - Ck dQ7j^y - + V • [(/x + cr k p t )Vk] (1) 

|p(kL) + V • (pU(kL)) - C 01 ® P k - C 02 pk 3 / 2 - (2) 

+v • [(/x + <r^,/Xt)V (kL)] 


The production of turbulent kinetic energy is stress-based, Eq. 3. and is limited in both the 
k and kL equations, Eq. 4. 


P — Tjj ^ ^ , Tjj — 2/i r Sij ^ pk(5y , Sy 

Pt = mi „(p,2ocy/^) 

with the turbulent eddy viscosity computed using Eq. 5. 


M t = 


i /4 p(kL) 

kVa 


The functions are: 


Ck = c// 4 , c 01 = Cl - C2 = C 3 

1 + C d i? r _ pv^d 
U l + £ 4 ’ * 20 /x 

<9 2 Ui d 2 Ui 

dx\ 3a; 2 


L„a- — K 


U' 

U" 


U' = v/2SySy , U" = 


1 f chM 

2 \dxj 



(3) 


(4) 


(5) 

( 6 ) 

(7) 

( 8 ) 
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The second derivative expression of the velocity can be written out as: 

' d 2 u d 2 u d 2 u\~ ( d 2 v d 2 v d 2 v\~ ( d 2 w d 2 w d 2 w\' 

1 1 -i- 1 1 -i- 1 1 


U" = 


dx 2 dy 2 dz 2 


dx 2 dy 2 dz 2 


\ dx 2 dy 2 dz 2 ) 


A limiter is applied on L„£, 

Lvfc,min A L-ufc A htifc,max; 




(kL) 

kCn 


Li;fc,max — Ci2«:d/p 


/ p = min 


( Pk(kL) 

\ C 3 / 4 pk 5 / 2 ’ 


,0.5 ,1.0 


( 9 ) 

(10) 

( 11 ) 


The boundary conditions for the two turbulence variables k and (kL) along solid walls and 
recommended farfield boundary conditions for most applications are: 


i— WJ 
<x oo 1 


kwall — (k^)wall — koo — 9 X 10 
where a represents the speed of sound. The constants are: 


(kL)^ = 1.5589 x 10 - 6 ^°° a °° 


(12) 


Ok 

= 1 . 0 , 

O'(kL) = 1.0 

K 

= 0 . 41 , 

C M = 0.09 

Cl 

= 1 - 2 , 

C 2 = 0 . 97 , C 3 = 0.13 

Cn 

= 10 . 0 , 

C 12 = 1 . 3 , C dl = 4.7 


2.2 Jet Corrections for Free Shear and Compressibility 

Subsonic and supersonic jet flows are quite difficult to predict with most RANS turbulence 
models. For subsonic jets, turbulence models do not predict the correct mixing rate and 
predict core lengths that are either too long or too short. The base k-kL-MEAH2015 model 
also predicts too short of core length. 

However, it turns out that the k-kL-MEAH2015 model predicts the correct mixing rate 
and shorter core length. By modifying the k-equation diffusion coefficient, the core length 
is improved, using Eqs. 13 and 14. 


ok = / 2 CTk! + (1 - / 2 ) cr k2 , CTki = 1-0, cr k2 = 0.5, (13) 

h = tank (r 2 ) , r = max (CAjL < 14 > 

This modification does not affect attached flow simulations and is active only in the wake 
or shear flow regions. This is similar to the approach used by Menter in the SST model to 
switch between k-w and k -e via the / 2 function. 

Most turbulence models fail to predict high-speed shear flow, as the mixing is much slower 
than subsonic flow. Compressibility correction is the approach used by most turbulence 
models to improve this deficiency. We propose to use an approach similar to Wilcox’s 
compressibility with cut-off Mach number to activate the compressibility for supersonic flow 
and not affect subsonic shear flow, as listed in Eqs. 15 and 16. 

C k = C 3 / 4 (1 + f c ) , C 02 = & + 2.5C 3/4 / c (15) 

fc = 1.5 (M t 2 - M 2 ) H [M 2 - M 2 ] , M t = M 0 = 0.19 (16) 

This jet corrected model, including both free shear and compressibility correction terms, is 
termed k-kL-MEAH2015+J. 


5 



3 Computational Methods 
3.1 CFL3D Code 

CFL3D [11] is a structured-grid upwind multi-zone CFD code that solves the generalized 
thin-layer or full Navier-Stokes equations. In the current study, the full viscous terms are 
used for all computations. CFL3D can use point-matched, patched, or overset grids and 
employs local time-step scaling, grid sequencing and multigrid to accelerate convergence 
to steady state. A time-accurate mode is available, and the code can employ low-Mach 
number preconditioning for accuracy in computing low-speed steady-state flows. CFL3D is 
a cell-centered finite-volume method. It uses third-order upwind-biased spatial differencing 
on the convective and pressure terms, and second-order differencing on the viscous terms; 
it is globally second-order accurate. Roe’s flux difference-splitting method [12] is used to 
obtain fluxes at the cell faces. The solution is advanced in time with an implicit approxi- 
mate factorization method. For each loosely coupled iteration, the mean flow equations are 
advanced in time with the eddy-viscosity fixed then the turbulence model is advanced in 
time with the mean flow solution fixed. A wide variety of turbulence models are available 
in the code, including linear eddy-viscosity and nonlinear models. 


3.2 Fun3D Code 

Fun3D is an unstructured three-dimensional, implicit, Navier-Stokes code. Roe’s flux dif- 
ference splitting [12] is used for the calculation of the explicit terms. Other available flux 
construction methods include HLLC [13], AUFS [14], and LDFSS [15]. The default method 
for calculation of the Jacobians is the flux function of van Leer [16], but the method by 
Roe and the HLLC, AUFS and LDFSS methods are also available. The use of flux limiters 
are mesh and flow dependent. Flux limiting options include MinMocl [17] and methods by 
Barth and Jespersen [18] and Venkatakrishnan [19]. Other details regarding Fun3D can be 
found in Anderson and Bonhaus [20] and Anderson et al. [21], as well as in the extensive 
bibliography that is accessible at the Fun3D Web site, http://fun3d.larc.nasa.gov. 


4 Test Case Descriptions 

The five geometries used in the verification and validation of the turbulence model are de- 
scribed in sub-sections 4.1 through 4.6. Table 1 lists relevant aspects of the simulation for 
each case. The flat plate and planar shear flows are low speed with little, if any, compress- 
ibility effects. The transonic bump flow has both flow separation and re-attachment points. 
The jet flows display both low and high compressibility characteristics. Comparisons are 
made with historic, canonical data or experimental results where available. 

Table 1. Test Cases. 


Geometry 

Grid 

Flow physics 

Flat plate 

2D 

wall bounded, attached, steady 

Planar shear 

2D 

free shear, steady 

Transonic bump 

Axisymmetric 

wall bounded, separated, unsteady 

Subsonic jet 

Axisymmetric 

free shear, low speed 

Supersonic jet 

Axisymmetric 

free shear, compressible 
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Figure 1. Flat plate geometry and boundary conditions. 


4.1 Zero Pressure Gradient Flat Plate 

Five successively refined grid levels are used to determine the numerical discretization error 
and the results from both codes are expected to converge to similar values. Different flow and 
turbulent quantities from each code are compared to additionally verify the implementation, 
along with available theoretical and experimental data. 

Figure 1 shows the sketch of the flat plate test case with boundary conditions used in 
this analysis. This is a subsonic, M m = 0.2 case at Re = 5 million per unit length. This is 
a verification and validation case. The following plot in Fig. 2 shows the convergence of the 
wall skin friction coefficient at x = 0.97 using five levels of grid size with the k-kL turbulence 
model. Each coarser grid is exactly every-other-point of the next finer grid, ranging from 
the super fine grid of ( 544 x 384 ) cells to the very coarse grid of ( 34 x 24 ) cells. In 
the plot, the x-axis is plotting yj (1/7V), which is proportional to grid spacing ( h ). At the 
left of the plot, h = 0 represents an infinitely fine grid. The difference in the skin friction 
coefficient between the coarsest and finest grid is less than 0.0002. Fun 3D and CFL3D 
converge toward the same result as the grid is refined. 

The skin friction coefficient, using the k-kL turbulence model on the finest grid, ( 544 
x 384 ) cells, over the entire plate, varies with respect to momentum thickness Reynolds 
number, as shown in Fig. 3. The results are computed using CFL3D and lie within the 
range of different correlation as compared with Karman-Schoenherr theory [22]. Using the 
finest grid results, Fig. 4 shows u + velocity with respect to y + predicted well compared with 
Cole’s theory [23], and also within the range of different correlations. Turbulence kinetic 
energy, length scale (kL) and viscosity are other quantities used to validate the quality 
of turbulence model predictions. Figure 5 shows comparisons of the predicted turbulence 
quantities results from Fun3D and CFL3D codes. Results from the two codes on this grid 
are indistinguishable, with the exception of a very small differences close to the wall for 
the variable kL. This could be the result of the difference in solver strategy. Fun3D is a 
node-centered code, whereas CFL3D is a cell-centered code. Figure 6 shows the comparisons 
between k-kL-MEAH2015 and SST-V turbulence models using CFL3D. Both models yield 
comparable results. 
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Figure 2. Skin friction coefficient with grid refinement, flat plate, x = 0.97, = 0.2, 

Re L = 5 million, k-kL-MEAH2015. 



4000 6000 8000 10000 12000 14000 

Re„ 


Figure 3. Skin friction correlation. — , Karman-Schoenherr [22]; k-kL-MEAH2015. 
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Figure 4. Velocity profile comparison of Cole’s theory with k-kL turbulence model, Re« = 
10,000, Coles [23], 


4.2 Planar Shear 

The planar shear case focuses on the development of the free shear layer following the passing 
of two different streams over a thin plate. The smaller inner stream has Mach number near 
M = 0.5, whereas the outer larger stream has a Mach number near M = 0.25. The Reynolds 
number is ReL = 50,000 based on unit length of the grid (L = 1). The computational 
domain extends from — 10 < x < 200, and 0 < y < 100. The separating plate extends from 
— 10 < x < 0 at y = 0.5. In terms of the plate, the reference length is 10 units. Both the 
lower and upper boundaries are taken to be symmetry planes. Figure 7 shows the close-up 
of the case with boundary conditions listed. 

This is a verification case. The plot in Fig. 8 shows the convergence of the drag coefficient 
on the divided plate using five levels of grid size with the k-kL-MEAH2015 turbulence model. 
Each coarser grid is exactly every-other-point of the next finer grid, ranging from the super 
fine grid of three blocks ( 128 x 256, 128 x 256, 512 x 512 ) containing 327,680 cells to the 
very coarse grid of ( 8 x 16, 8 x 16, 32 x 32) containing 1280 cells. In the plot, the x-axis 
plots yj (1/iV), which is proportional to grid spacing (h). The difference in the total skin 
friction coefficient comparing the coarsest and finest grid is less than 0.0002. 

Fun 3D and CFL3D converge toward similar values as the grid is refined. Figure 9 shows 
comparisons of the predicted turbulence quantities results from Fun 3D and CFL3D codes. 
Results from the two codes on the finest grid are essentially indistinguishable for k, kL and 
turbulence viscosity. The fully developed turbulent shear flow exhibits self-similar behavior 
downstream. This can be achieved by normalizing the velocity and normal distance. The 
velocity can be normalized as ( u — u e ) / ( u m — u e ), where u e is the velocity at the edge of the 
outer stream, and it m is the peak (centerline) velocity. When plotted against y/b, where b 
is the half width (location where ( u — u e ) is half of (u m — u e ), the results can be compared 
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k/a* 



(a) Turbulent kinetic energy. 


(b) kL. 



M-/H, 


(c) Turbulent eddy viscosity. 

Figure 5. Code-to-code comparison, flat plate, = 0.2, ReL= 5 million, x = 0.95, k-kL- 
MEAH2015. 
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Figure 6. Comparison of velocity profiles, flat plate, = 0.2, ReL= 5 million, x = 0.97. 


2 - 


1.5 - 


0.5 


- 0.5 


p/p„ = 1 .046 

VC = i -0 

1 quantity from interior 


adiabatic solid wall 

/ 


p/p„= 1.1862 V 

T/C = 1 -05 adiabatic solid wail 

1 quantity from interior 


-10 


t 

symmetry 

i I i 

o 


X 


Figure 7. Sketch of geometry and boundary conditions. 
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Figure 8. Effect of grid refinement on drag coefficient, planar shear, M™ = 0.2, Rer = 
50,000, k-kL-MEAH2015. 


to the experimental data of Bradbury and Riley [24]. In Fig. 10a, CFL3D results are taken 
from the three x-locations x = 29.2468, x = 64.2188, and x = 95.501. The CFD results 
using the k-kL-MEAH2015 turbulence model are approximately self-similar, very tight, and 
agree well with the experiment. Similar comparisons are shown in Fig. 10b using SST-V, 
which did not agree as well with the data as k-kL-MEAH2015. 

4.3 Axisymmetric Transonic Bump 

Figure 11 shows the sketch of the axisymmetric transonic bump case with boundary con- 
ditions used in this analysis. This is a transonic, = 0.875 case at ReL = 2.763 million 
based on L = chord length. The purpose here is to provide a validation case that establishes 
the model’s ability to predict separated flow. For this particular axisymmetric transonic 
bump case, the experimental data are from Bachalo and Johnson [25]. The experiment 
utilized a cylinder of 0.152 m diameter in a closed return, variable density, and continuous- 
running tunnel with 21% open porous-slotted upper and lower walls. The boundary layer 
incident on the bump was approximately 1 cm thick. The bump chord was 0.2032 m. In 
the experimental case, with a freestream Mach number of 0.875, the shock and trailing-edge 
adverse pressure gradient results in flow separation with subsequent reattached downstream 
of the bump. Figure 13 shows surface pressure and velocity as a function of x, comparing 
experimental data and results from Fun3D and CFL3D using the k-kL-MEAH2015 turbu- 
lence model. Both CFD codes produced similar results. In general, surface pressure is in 
good agreement with data. The velocity distributions also agree well with experimental data 
for all the stations, and are generally better than the results generated using SST, shown in 
Fig. 12, especially at the x/c = 0.668 location. The majority of RANS turbulence models fail 
to predict the correct reattachment of this separated flow case, with the separation bubble 
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Figure 9. Flow profiles, planar shear, x = 29.2468, = 0.2, Re L = 50,000, k-kL- 

MEAH2015. 
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Figure 10. Self-similar velocity profiles, planar shear, = 0.2, Re L = 50,000. symbols, 
Data-Bradbury & Riley [24]; lines, u m /u e ~ 0.5, CFL3D, k-kL-MEAH2015. 
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Figure 11. Sketch of transonic bump configuration. 
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x/c 

Figure 12. Comparison of surface static pressure coefficients, = 0.875. symbols, Data- 
Bachalo & Johnson [25]; lines, k-kl-MEAH2015. 


size over predicted by at least 25%. 

The prediction of flow separation and reattachment location using k-kL-MEAH2015 and 
SST is tabulated in Table 2. The bubble size predicted with k-kL-MEAH2015 is reduced to 
less than 8% compared with 27% from the SST simulation. 

Table 2. Transonic bump. 


Data 

Experiment 

k-kL 

SST 

Separation location 

0.70 

0.68 

0.65 

Reattachment location 

1.10 

1.11 

1.16 

Bubble size error (%) 

- 

7.50 

27.5 
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u/U, 


(a) k-kL-MEAH2015 turbulence model. 



u/U. 


(b) SST turbulence model. 


Figure 13. Velocity profile comparisons at = 0.875 between Bachalo & Johnson [25] 
(symbols) and CFD (lines). 


4.4 Subsonic Jet: Bridges 

Axisymmetric subsonic and near sonic jet cases are used to validate the prediction of jet 
flows. The experiment measured velocities as well as turbulence quantities downstream of 
the jet exit using Particle Image Velocimetry (PIV) [26]. Velocity profiles of interest are 
measured at the centerline ( y = 0 ) as well as other x-locations. In the present report, we 
are comparing the turbulence model results with the centerline values for velocity data. The 
grid used was obtained from [9] and was converted to unstructured hex grid with 145,561 
nodes, and solved only using the Fun 3D code. In the experiment, the axisymmetric jet 
exits into quiescent (non-moving) air at two nozzle exit Mach conditions, M cx i tj acoustic = 
ujet / aoo = 0.51 and 0.9. However, because simulating flow into totally quiescent air is 
difficult to achieve for some CFD codes, here the solution is computed with a very low 
background ambient condition of = 0.01, moving left-to-right in the same direction as 
the jet). Figure 14 shows the grid and flow setup for the subsonic jet case. The ability of 
the k-kL-MEAH2015-J turbulence model in predicting jet flow, compared with the basic 
k-kL-MEAH2015 and SST turbulence models, is shown in Figs. 15 and 16. 


[h] 

Table 3. Subsonic jet conditions. 


Set Point 

Miich-exit , acoustic 

NPR 

NTR 

Tjet, static/ Tqq 

3 

0.51 

1.197 

1.0 

0.950 

7 

0.90 

1.861 

1.0 

0.835 


For Set Point 3 in Fig. 15, the jet core length and rate of decay are better predicted 
when using the k-kL-MEAH2015+J model, as shown in Fig. 15a. In particular, the jet core 
length is in better agreement with experiment than the much longer core predicted using 
the SST turbulence model. Very good comparisons of k-kL-MEAH2015+J are shown in 
Fig. 15b for the velocity variations with radial direction at different x locations. Figure 16 
shows results for Set Point 7. Again the k-kL-MEAH2015+J turbulence model predicts the 
jet flow better than the basic k-kL-MEAH2015 and SST turbulence models. As shown in 
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Figure 14. ARN1 geometry and boundary conditions. 


Fig. 16a, the jet core length and rate of decay is in better agreement with experiment. Very 
good comparisons of k-kL-MEAH2015+J are shown in Fig. 16b for the velocity variations 
with radial direction at different x locations. 

4.5 Supersonic Jet: Seiner 

In the subsonic jet cases above, the first part of the jet correction terms (free shear correction) 
improved the prediction of the k-kL-MEAH2015 turbulence model by shifting the jet core 
downstream. The results are in generally good agreement with experimental data and are 
better than SST predictions. The second part of the jet correction terms (compressibility 
correction) has no effect in either of the subsonic jet cases, because the turbulence Mach 
number is smaller than the cutoff at 0.19. The term will be activated around jet Mach 
number of 1.5. 

The supersonic jet cases, Seiner [27] and Eggers [28], are frequently used to validate 
turbulence models for high-speed flow. RANS models are generally not capable of predicting 
the mixing rate and core length of high-speed jet flow without adding a compressibility 
correction. Both cases are for nozzles operating near design condition and fully expanded. 
However, because flow into quiescent air is difficult to achieve for some CFD codes for such 
high shear flow cases, the simulations are run with low background ambient conditions of 
M^ = 0.05. Three grid levels are utilized to verify the implementation of the proposed 
modification and show grid convergence. The first case is Seiner’s Mach 2.0 C-D nozzle 
operating at NPR = 7.824. The grid was obtained from the NPARC website [29] and 
converted to an unstructured hex grid. The grid was modified to provide 3 grid levels 28433, 
112865 and 449729 nodes. The medium grid is shown in Fig. 17. These cases were run only 
using the Fun3D code. Figure 18 shows the grid effects on the prediction of centerline 
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Data, Bridges 


- k-kL-MEAH2015 


k-kL-MEAH2015+J FUN3D 


- SST 



(a) Streamwise centerline velocity. 



Figure 15. Comparison of jet velocity data, ARN1, Set Point 3. symbols, Data-Bridges [26]; 
lines, FUN3D. 


Mach number and total pressure using k-kL-MEAH2015+J turbulence model. All three- 
grid levels give similar results and the difference is very small between the medium and fine 
grid levels. As shown in Fig. 19 using medium grids, both the basic k-kL-MEAH2015 and 
SST turbulence models underpredict the jet core length for centerline Mach number and 
total pressure. This behavior is observed when using many other RANS turbulence models. 
The k-kL-MEAH2015+J model improves the results dramatically and achieves very good 
agreement with experimental data. 

4.6 Supersonic Jet: Eggers 

The second case is Eggers’ supersonic C-D nozzle test [28] and is sketched in Fig. 20. The 
grid was obtained from [30] and converted to unstructured hex grid. The nozzle geometry 
produces an exit Mach number of 2.2 when fully expanded and was set to the operating 
condition of NPR = 11.03. The grid was modified to produce 3 grid levels consisting of 
21432, 84563 and 335925, nodes each and all simulations were run with Fun3D. Figure 21 
shows the grid effects on the prediction of centerline Mach number velocity variation at 
x/r=28.93 using k-kL-MEAH2015+J turbulence model. All three-grid levels give similar 
results and the difference is very small between the medium and fine grid levels. Using the 
medium grid, both the basic k-kL-MEAH2015 and SST turbulence models underpredict the 
jet core length as shown in Fig. 22 for centerline velocity. As seen in the previous section, the 
jet correction model with the baseline k-kL-MEAH2015 improves the comparison with the 
experimental results. Similarly, the k-kL-MEAH2015+J model produces better matching 
for the variations of axial velocity with radial direction at different x locations as shown in 
Fig. 23. 
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n/n 




(a) Streamwise centerline velocity. 


(b) Radial velocity profiles, k-kL-MEAH2015+J. 


Figure 16. Comparison of jet velocity data, ARN1, Set Point 7. symbols, Data-Bridges [26]; 
lines, FUN3D. 
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Figure 17. Axisymmetric nozzle geometry and conditions, Seiner [27]. 
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Figure 18. Effect of grid refinement, NPR = 7.824, Mach 2. jet. symbols, Data-Seiner [27]; 
lines, FUN3D, k-kL-MEAH2015+J. 



Figure 19. Effect of jet correction on flowfield, NPR = 7.824, Mach 2. jet. symbols, 
Data-Seiner [27]; lines, FUN3D. 
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Figure 20. Geometry and boundary conditions, Eggers [28]. 




(a) Streamwise centerline velocity. 


(b) Radial velocity profile. 


Figure 21. Effect of grid refinement on velocities, Mach 2. jet. symbols, Data-Eggers [28]; 
lines, FUN3D, k-kL-MEAH2015+J. 
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Figure 22. Effect of jet correction on centerline velocity, Mach 2. jet. symbols, Data- 
Eggers [28]; lines, FUN3D. 
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Figure 23. Effect of jet correction on velocity profiles, Mach 2. jet. symbols, Data- 
Eggers [28]; lines, FUN3D. 
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5 Summary 


With the exception of Reynolds stress, the turbulent kinetic energy is exclusively used as 
one of the transport equations in multi-equation turbulence models. The foundation of 
this equation is well established and accepted. The scale-determining equation, though, is 
considered the weakest link, even when full Reynolds stress and hybrid RANS/LES for- 
mulations are considered. The first objective of this report is to validate and verify the 
implementation of the newly developed turbulence model k-kL-MEAH2015 in both CFL3D 
and Fun3D CFD codes. Rotta shows that a reliable scale-determining equation can be 
formed in a transport equation for the turbulent length scale, L. Rotta’s equation is well 
suited for a term-by-term modeling and shows some interesting features compared to other 
approaches. The most important difference is that the formulation leads to a natural inclu- 
sion of higher order velocity derivatives into the source terms of the scale equation. Five 
steady-state test cases were computed and presented, including an attached flow flat plate, 
2D bump, separated flow, shear flow, and jet flow cases. 

The flat plate and 2D bump are verification cases that use 5 grid levels in a conver- 
gence study. All the verification cases show grid convergence and independent solutions. 
The implementation of k-kL-MEAH2015 in both CFL3D and Fun3D is verified from these 
results. The flat plate case is also used as a validation case and compared with theoreti- 
cal data. The results from these cases are compared with theory and experimental data, 
as well as with results using the SST turbulence model. They demonstrate that the k-kL- 
MEAH2015 model has the ability to produce results similar or better than SST results. The 
size of the separation bubble (separation and re-attachment locations) is better predicted 
by k-kL-MEAH2015 model for the separated axisymmetric transonic bump case. 

Additionally, the k-kL-MEAH2015 model with jet correction gives better agreement 
with experimental data than the basic k-kL-MEAH2015 and SST models for the subsonic, 
near-sonic and supersonic jet cases. 


References 

1. Launder, B. E. and Spalding, D. B.: Lectures in mathematical models of turbidence, 
Academic Press, London, 1972. 

2. Rotta, J.: “Statistische Theorie nichthomogener Turbulenz,” Zeitschrift fiir Physik, Vol. 
129, No. 6, 1951, pp. 547 -572. 

3. Menter, F. R.; Egorov, Y. and Rusch, D.: “Steady and Unsteady Flow Modeling Using 
the k-\/kL model,” Turbulence, Heat and Mass Transfer, 5th International symposium, 
Proceedings , eds. K. Hanjalic, Y. Nagano and S. Jakirlic, Dubrovnik, Croatia, 2006, pp. 
403-406. 

4. Menter, F. and Egorov, Y.: “The Scale- Adaptive Simulation Method for Unsteady 
Turbulent Flow Predictions. Part 1: Theory and Model Description,” Flow, Turbulence 
and Combustion , Vol. 85, No. 1, 2010, pp. 113-138. 

5. Egorov, Y.; Menter, F.; Lechner, R. and Cokljat, D.: “The Scale-Adaptive Simulation 
Method for Unsteady Turbulent Flow Predictions. Part 2: Application to Complex 
Flows,” Flow, Turbulence and Combustion, Vol. 85, No. 1, 2010, pp. 139-165. 

6. Menter, F. R.: “Improved Two-Equation k -oj Turbulence Models for Aerodynamic 
Flows,” NASA TM-103975, Oct. 1992. 


24 



7. Abdol-Hamid, K. S.: “Assessments of k-kl Turbulence Model Based on Menter’s Modi- 
fication to Rotta’s Two-Equation Model,” International Journal of Aerospace Engineer- 
ing, Vol. 2015, 2015, pp. 1-18. 

8. Abdol-Hamid, K. S.; Pao, S. P.; Hunter, C.; Deere, K.; Massey, S. and Elmiligui, A.: 
“PAB3D: It’s History in the Use of Turbulence Models in the Simulation of Jet and 
Nozzle Flows,” AIAA Paper 2006-0489, Jan. 2006. 

9. http://turbmodels.larc.nasa.gov, accessed: 2015-09-01. 

10. Rumsey, C.; Smith, B. and Huang, G.: “Description of a Website Resource for Turbu- 
lence Model Verification and Validation,” AIAA Paper 2010-4742, June 2010. 

11. Krist, S. L.; Biedron, R. T. and Rumsey, C. L.: “CFL3D User’s Manual (Version 5.0),” 
NASA TM-1998-208444, June 1998. 

12. Roe, P. L.: “Approximate Riemann Solvers, Parameter Vectors, and Difference 

Schemes,” J. Comp. Phys., Vol. 43, 1981, pp. 357-372. 

13. Batten, P.; Clarke, N.; Lambert, C. and Causon, D.: “On the Choice of Wavespeeds for 
the HLLC Riemann Solver,” SIAM J. Sci. Comput., Vol. 18, 1997, pp. 1553-1570. 

14. Sun, M. and Takayama, K.: “Artificially Upwind Flux Vector Splitting Scheme for the 
Euler Equations,” J. Comp. Phys., Vol. 189, 2003, pp. 305-329. 

15. Edwards, J.: “A Low-Diffusion Flux-Splitting Scheme for Navier Stokes Calculations,” 
AIAA Paper 1996-1704, May 1996. 

16. van Leer, B.: “Towards the Ultimate Conservative Difference Schemes V. A second 
order sequel to Godunov’s Method,” J. Comp. Phys., Vol. 32, 1979, pp. 101 -136. 

17. Roe, P. L.: “Characteristic-Based Schemes for the Euler Equations,” Annual Review of 
Fluid Mechanics, Vol. 18, 1986, pp. 337-365. 

18. Barth, T. and Jespersen, D.: “The Design and Application of Upwind Schemes on 
Unstructured Meshes,” AIAA Paper 1989-0366, Jan. 1989. 

19. Venkatakrishnan, V.: “Convergence to Steady State Solutions of the Euler Equations 
on Unstructured Grids with Limiters,” J. Comp. Phys., Vol. 118, 1995, pp. 120-130. 

20. Anderson, W. and Bonhaus, D.: “An Implicit Upwind Algorithm for Computing Tur- 
bulent Flows on Unstructured Grids,” Computers and Fluids, Vol. 23, No. 1, 1994, pp. 
1 - 22 . 


21. Anderson, W.; Rausch, R. and Bonhaus, D. L.: “Implicit/Multigrid Algorithms for 
Incompressible Turbulent Flows on Unstructured Grids,” J. Comp. Phys, Vol. 128, 
1996, pp. 391-408. 

22. Schoenherr, K. E.: Resistance of Flat Surfaces Moving Through a Fluid, Ph.D. thesis, 
Johns Hopkins University, 1932. 

23. Coles, D.: “The law of the wake in the turbulent boundary layer,” J. Fluid Mech., 
Vol. 1, 1956, pp. 191-226. 

24. Bradbury, L. J. S. and Riley, J.: “The spread of a turbulent plane jet issuing into a 
parallel moving airstream,” J. Fluid Mech., Vol. 27, 1967, pp. 381-394. 


25 



25. Bachalo, W. D. and Johnson, D. A.: “Transonic, turbulent boundary-layer separation 
generated on an axisynnnetric flow model,” AIAA Journal, Vol. 24, No. 3, 1986, pp. 
437-443. 

26. Bridges, J. and Brown, C.: “Parametic Testing of Chevrons on Single Flow Hot Jets,” 
AIAA Paper 2004-2824, 2004. 

27. Seiner, J. M.; Ponton, M. K.; Jansen, B. J. and Lagen, N. T.: “The Effects of Tempera- 
ture on Supersonic Jet Noise Emission,” DGLR / AIAA Aeroacoustics Conference, 14th, 
Proceedings, Vol. 1, (A93-19126 05-71), Aachen, GERMANY, May 1992, pp. 295-307. 

28. Eggers, J.: “Velocity Profiles and Eddy Viscosity Distributions Downstream of a Mach 
2.22 Nozzle Exhausting to Quiescent Air,” NASA TN D-3601, September 1966. 

29. http://www.grc.nasa.gov/WWW/wind/valid/seiner/Seiner .html, accessed: 2015- 

09-10. 

30. http://www.grc.nasa.gov/WWW/wind/valid/axinoz/axinoz.html, accessed: 2015- 

09-10. 


26 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 
Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 
shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 


1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

01-11 - 2015 Technical Memorandum January 2005 - December 2008 

4. TITLE AND SUBTITLE 

Verification and Validation of the k-kL Turbulence Model in FUN3D and 
CFL3D Codes 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Abdol-Hamid, Khaled S.; Carlson, Jan-Renee; Rumsey, Christopher L. 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

109492.02.07.01.01.01 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

NASA Langley Research Center 
Hampton, VA 23681-2199 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

L-20609 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSOR/MONITOR'S ACRONYM(S) 

NASA 

11. SPONSOR/MONITOR'S REPORT 
NUMBER(S) 

NASA-TM-20 15-2 18968 


12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category 02 

Availability: NASA STI Program (757) 864-9658 

13. SUPPLEMENTARY NOTES 


14. ABSTRACT 

The implementation of the k-kL turbulence model using multiple computational fluid dynamics (CFD) codes is reported herein. The k-kL model is a 
two-equation turbulence model based on Abdol-Hamid's closure and Menter's modification to Rotta's two-equation model. Rotta shows that a reliable transport 
equation can be formed from the turbulent length scale L, and the turbulent kinetic energy k. Rotta's equation is well suited for term-by-term modeling and 
displays useful features compared to other two-equation models. An important difference is that this formulation leads to the inclusion of higher order velocity 
derivatives in the source terms of the scale equations. This can enhance the ability of the Reynolds-averaged Navier-Stokes (RANS) solvers to simulate 
unsteady flows. The present report documents the formulation of the model as implemented in the CFD codes FUN3D and CFL3D. Methodology, verification 
and validation examples are shown. Attached and separated flow cases are documented and compared with experimental data. The results show generally very 
good comparisons with canonical and experimental data, as well as matching results code-to-code. The results from this formulation are similar or better than 
results using the SST turbulence model. 


15. SUBJECT TERMS 

CFD; Navier-Stokes; Reynolds averaged; Turbulence models; Validation; Verification 


16. SECURITY CLASSIFICATION OF: 


a. REPORT b. ABSTRACT c. THIS PAGE 


17. LIMITATION OF 
ABSTRACT 


18. NUMBER 
OF 

PAGES 


19a. NAME OF RESPONSIBLE PERSON 

STI Help Desk (email: help@sti.nasa.gov) 

19b. TELEPHONE NUMBER (Include area code) 


(757) 864-9658 


Standard Form 298 (Rev. 8-98) 

Prescribed by ANSI Std. Z39.18 


u 


u 


u 


uu 


31 



































